#include<iostream>
#include"mpi.h"
#include<immintrin.h>
#include<sys/time.h>
#include<cmath>

using namespace std;

int myid, numprocs, namelen, cpuname, masterNode;//进程标识号与进程总数，可在任意函数中访问
MPI_Status status;
char processor_name[MPI_MAX_PROCESSOR_NAME];
int partRow, up, down;//各进程分配到的行数，相邻的上下两个进程
const int sca = 16384;//问题规模

void Lu(int N, float** A);//普通消元算法
void printm(int N, float** A);
void generatem(int N, float** A);

void LU_mpi_cir(int argc, char* argv[], float** A) {  //等步长的循环划分
	double start_time = 0;
	double end_time = 0;
	//MPI_Init(&argc, &argv);
	int total = 0;
	int rank = 0;

	MPI_Status status;
	MPI_Comm_size(MPI_COMM_WORLD, &total);
	MPI_Comm_rank(MPI_COMM_WORLD, &rank);
	if (rank == 0) {  //0号进程初始化矩阵//generatem(sca, A);
		//printm(sca, A);
		for (int j = 1; j < total; j++) {
			for (int i = j; i < sca; i += total) {
				MPI_Send(&A[i][0], sca, MPI_FLOAT, j, 1, MPI_COMM_WORLD);//1是初始矩阵信息，向每个进程发送数据
			}
		}
	}
	else {
		for (int i = rank; i < sca; i += total) {
			MPI_Recv(&A[i][0], sca, MPI_FLOAT, 0, 1, MPI_COMM_WORLD, &status);
		}
	}
	MPI_Barrier(MPI_COMM_WORLD);
	//start_time = MPI_Wtime();
	for (int k = 0; k < sca; k++) {
		if (k % total == rank) {
			for (int j = k + 1; j < sca; j++) {
				A[k][j] = A[k][j] / A[k][k];
			}
			A[k][k] = 1.0;
			for (int j = 0; j < total; j++) { //
				if (j != rank)
					MPI_Send(&A[k][0], sca, MPI_FLOAT, j, 0, MPI_COMM_WORLD);//0号消息表示除法完毕
			}
		}
		else {
			int src = k % total;

			MPI_Recv(&A[k][0], sca, MPI_FLOAT, src, 0, MPI_COMM_WORLD, &status);
		}
		int begin = k;
		while (begin % total != rank)
			begin++;
		for (int j = begin; j < sca; j += total) {
			for (int i = k + 1; i < sca; i++) {
				A[i][j] = A[i][j] - A[i][k] * A[k][j];
			}
		}
	}
	MPI_Barrier(MPI_COMM_WORLD);	//各进程同步
	MPI_Finalize();
}


int main(int argc, char** argv)
{
	struct timeval sign1, sign2, sign3, sign4;

	const int siz = sca;
	float** A = new float* [siz];
	generatem(siz, A);
	//printm(siz, A); 
	
	gettimeofday(&sign1, NULL);
	LU_mpi_cir(argc, argv, A);
	
	gettimeofday(&sign2, NULL);
	Lu(siz, A);
	
	gettimeofday(&sign3, NULL);

	
	cout << "m+s算法时间/us： " 
		<< (sign2.tv_sec*1000*1000+sign2.tv_usec)
			-(sign1.tv_sec*1000*1000+sign1.tv_usec);
	cout << "m算法时间/us： " 
		<< (sign3.tv_sec*1000*1000+sign3.tv_usec)
			-(sign2.tv_sec*1000*1000+sign2.tv_usec);
	cout << "\n";

	return 0;
}

void Lu(int N, float** A)//普通消元算法
{
	for (int k = 0; k < N; k++)
	{
		for (int j = k + 1; j < N; j++)
		{
			A[k][j] = A[k][j] / A[k][k];
		}//归一化
		A[k][k] = 1.0;
		for (int i = k + 1; i < N; i++)
		{
			for (int j = k + 1; j < N; j++)
			{
				A[i][j] = A[i][j] - A[i][k] * A[k][j];
			}
			A[i][k] = 0;
			//cout << i << ":\n";
			//print(N, A);
		}
	}
}
void printm(int N, float** A)
{
	for (int i = 0; i < N; i += 1)
	{
		for (int j = 0; j < N; j += 1)
		{
			cout << A[i][j] << ' ';
		}
		cout << '\n';
	}
}
void generatem(int N, float** A)
{
	//A = new float* [N];
	for (int i = 0; i < N; i += 1)
	{
		A[i] = new float[N];
		for (int j = 0; j < N; j += 1)
		{
			//cin >> A[i][j];
			A[i][j] = (i + 1) * (j + 1);
			if (i > j)
				A[i][j] = 0;
		}
	}
	//printm(20, A);//只打印20*20看看效果
	for (int i = 0; i < N; i += 1)
	{
		for (int j = 0; j < N; j += 1)
		{
			A[i][j] += A[i / 2][j];
		}
	}//初始化乱排
}

